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Abstract 

Fractional generalizations of the Poisson process and branching Furry process are con- 
sidered. The link between characteristics of the processes, fractional differential equations 
and Levy stable densities are discussed and used for construction of the Monte Carlo 
algorithm for simulation of random waiting times in fractional processes. Numerical cal- 
culations are performed and limit distributions of the normalized variable Z = N/ (N) are 
found for both processes. 
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1 Introduction 

The Poisson process is the simplest but the most important model for physical and other appli- 
cations. Its main properties (absence of memory and jump-shaped increments) model a large 
number of natural and social processes. Basic equations of theoretical physics (Schroedinger's, 
Pauly's, Dirac's and other equations) are derived in frame of axioms of the Poisson process. 
These equations describe fundamental processes on a microscopic physical level. 

When investigating complex macroscopic systems, we can observe another kind of behavior 
showing the presence of memory. 

There exist a few fractional generalizations of the Poisson process [Repin & Saichev, 2000, 
Jumarie Guy, 2001, Wang Xiao-Tian & Wen Zhi-Xiong, 2003, Wang Xiao-Tian et al, 2006, 
Laskin, 2003]. We consider here the fractional Poisson process (fPp), introduced by the waiting 
time distribution density ip v {t) with the Laplace transform 

oo 

&(A)= / 'e- xt Mt)dt=^-- (1) 
o n 

The density is characterized by fractional exponent v G (0, 1], being the order of the fractional 
differential equation describing this process. When v = 1, the fPp becomes the standard 
Poisson process, 

&(A) = -^ T , Mt)=^; 
H + A 

when v < 1, the fPp displays the property of memory, namely, the correlations between events 
in non-overlapping time intervals arise. Thus, the memory is a function of the parameter v that 
makes the fPp very attractive for investigating the mathematical phenomenology of memory. 

We consider in this work the main properties of the fPp, reformulate them in terms of 
alpha-stable densities , construct an algorithm for simulation of interarrival time, apply it to 
Monte Carlo simulation of the fPp, and find limit distributions. On the base of this approach 
we will introduce also the fractional generalization of the simplest branching process and find 
the correspondent distributions. But first of all we'll bring some physical reasons of our interest 
to the fPp. 



2 Physics of fPp 



For the sake of clearness, let us talk about number of events N(t) as about coordinate of a 
particle performing spasmodic motion in a given direction. The particle, appeared at the origin 
at time t = 0, stays there some random time Ti, then makes a jump to the position N — 1, 
stays there random time T 2 , then the second jump to N = 2 and so forth. Suppose that all 
Tj are independent and identically distributed random variables. Let Q(t) = P(T > t). If 
Q(t) = e _/1 *, fj, > 0, then N(t) is the Poisson process with rate fi. If the distribution of T is not 
an exponential one, but (T) < oo, the process is not a Poisson one, but at large scales it looks 
like a Poisson process and could be called the asymptotically Poisson process. As in the first 
case, the motion of the particle considered at large scales will seem to be almost regular and 
uniform with the "velocity" 1/ (T). There are no special asymptotical properties which appear 
here. 

The situation changes when Q(t) has a power law long tail 

Q(t) oc r u , < v < 1, (2) 

and an infinite mathematical expectation (T). Namely this type of traps leads to a fractional 
generalization of the Poisson process. 

The physical grounds of (j2J) has been discussed in a number of works. The first of inter- 
pretation of (J2J) has been done by Tunaley [1972] on the base of the following simple jump 
mechanism. The process goes in an insulator containing randomly distributed point (of small 
size) traps with exponentially distributed waiting times: P{T > t\9} = exp(— tjff). Their mean 
value 9 is finite and linked with the random distant 5 to the nearest site in the direction of the 
applied field as follows [Harper, 1967]: 

9 = /3[exp( 7 5) - 1]. 

Here, 7 is a positive constant and (3 is inversely proportional to the applied potential gradi- 
ent, both are independent of the temperature of the sample. Taking for 5 the exponential 
distribution with the mean d, 

P{5 >x } = e~ x/d } 
we obtain the probability density for 9 in the following form 

P{9 >t} = P{5 > (I/7) ln(l + t//3)} = exp[-(l/ 7 d) ln(l + t/p)] = (1 + t/f3)- l '^ d \ 

Averaging these distribution over d 

00 

P{T > t} = - J P{T > t\t'}dP{9 > t'} ~ vT(u)(t/P)- v , t->oo, u= l/( 7 d) 


yields 

In [Scher & Montroll, 1975], it has been indicated that the dispersive behavior can be 
caused also by multiple trapping in a distribution of localized states. On the assumption that 
the localized states below the mobility edge fall off exponentially with energy, one can arrive 
at Eq. (j2J) with the exponent v depending on the temperature B = kT. In the frame of this 
model, called random activated energy model it is assumed, that 
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(i) the jump rate of a particle hopping over an energy barrier AE has the usual quasiclassical 
form 

W = Ae~ AE / e ; 

(ii) the conditional waiting time distribution corresponding to a given activation energy 
AE = e is exponential 

P{T > t\e} = e~ w{£)t ; 

and 

(in) the activation energy is a random variable with the Boltzman distribution density 

p(e) = (Gc)-^-/®-. 
Averaging over the activation energy results: 

oo oo 

P{T >t} = J?{T> t\e}p(e)de = J exp[-(Ae' £ ^ e )t] d(e' £ ^) = uT(u)(Aty v 
o o 

with v = @/@ c - Here, O c is the characteristic temperature defining the conduction band 
tail. For 6 < 6 C , thermalization dominates and the photoinjected carriers sink progressively 
deeper in increasing time; the transport becomes dispersive. For > C , the carriers remain 
concentrated near the mobility edge and the charge transit exhibit non- dispersive behavior. 
Consequently, the physical meaning of v is that it is representative of disorder: the smaller its 
value the more dispersive the transport. 

The fluorescence emission of single semiconductor colloidal nanocrystals such as CdSe/ZnS 
core-shell quantum dots (QDs), exhibits unusual intermittency behavior [Shimizu et al. 2001]. 
Under laser illumination, nanocrystals blink: QDs jumps between bright and dark states. Ex- 
perimental investigations of blinking QD fluorescence showed that on- and off-periods are dis- 
tributed according to inverse power laws with exponents and respectively. The mechanism of 
QD's blinking is not yet quite understood. Thermally activated ionization, electron tunnelling 
through fluctuating barriers and some other mechanisms were suggested as alternative ones 
yielding power on- and off-time distributions. An emission wavelength of QD fluorescence is 
simply tuned by changing the size of the nanocrystal. This makes them promising as the active 
medium in light-emitting diodes or lasers. Fluctuations in the intensity of QD fluorescence can 
limit such applications. From theoretical results one follows that fluctuations in the case of 
power on- and off-intervals distributions stay constant and not decrease with time. 

These theoretical results being in accordance with numerous experimental data represent 
additional physical reasons for high attention to the fractional Poisson process, fractional dif- 
ferential equations and long memory phenomena models. 



3 Waiting time density 

The fPp waiting time density ip v (t) can be represented in three equivalent forms. The first of 
them, given in [Repin & Saichev, 2000], is 

1 * 



ipu(t) = -Je x (f) l/ (/2t/x)dx, 
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where 

, ( t \ = sin(^7r) 

{) 7r[e + ^- i/ + 2cos(z/7r)]' 
This form allows to find asymptotical expressions for small and large time, 

and to perform numerical calculations of the density (see Fig 1.). 




Figure 1: The fPp waiting time distribution densities ([3]) for fi = 1 and v = 0.1(0.1)1. 
The second form, obtained in [Laskin 2003], 

P(T >t) = E v (-nt v ), 

= ^ = ^~ l E„A-^) (3) 



uses the Mittag-Leffler functions 



z n 



^ T(an + f3) ' 

We present here the third form, which will serve as a basis for Monte Carlo simulation of 
fPp's. 

The complement cumulative distribution function P(T > t) can be represented in the form 

oo 

P(T >t) = J e-^ /T "g iu \T)dr } (4) 
o 

where g^ u \r) is the one-sided a-stable density (see for details [Uchaikin & Zolotarev 1999], 
[Uchaikin 2003], [Dubkov & Spagnolo 2005] and [Dubkov & Spagnolo 2007]). 
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Indeed, expanding the exponential function in (j2J) 



DC 



fc=0 



and making use of the formula for negative order moments of the a-st able-density 



CO 



T(l + ukY 

u 

we obtain 

?V > - E - E j^y - *<-*">■ 



k=0 n ^=0 



4 Simulation of waiting times 



The following result solves the problem of simulation of random waiting times. 
The random variable T determined above has the same distribution as 

I \nU\ 1/u 

where S{y) is a random variable distributed according to g^ u '(r) and U is independent of S(v), 
is a uniformly distributed in [0, 1] random variable. 

Making use of the formula of total probability, let us represent (J3J) in the following form 

oo 

P(T >t)= J P(T > t\T)g M (T)dT, (5) 

where 



o 



P(T > t|r) = e-^/^ (6) 
is the conditional distribution. This means that 

P(T > t\r) = P(U < e-^lT") = P > ^ , 

or 

= ^Ia— (?) 
Because r is a fixed possible value of S(u), we obtain for unconditional interarrival time 

d \lnU\^ 



The random variable 



£ I lnfAl 1 ^ sin(^ 2 )[sin((l - z/)^)] 1 /^ 1 

//V" [sin (7rf/ 2 )] V- [In f/ 3 ] ' [ ' 

where U\, U2 and U3 are independent uniformly distributed on [0,1] random numbers. This 
conclusion follows from the Kanter algorithm for simulating S(v) [Kanter, 1975]. 

Note that when v — > 1 this algorithm reduces to standard rule of simulating random numbers 
with exponential distribution: 



rp d_ 



\nU\ 
H 
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5 The nth arrival time distribution 

Let T^ n \ n = 1, 2, 3, . . . , be the nth arrival time of fPp 

T n = T (1) + T (2) + • • • + T (n) (9) 
and ip* n (t) be its probability density : 

n times 

Here, T^'s are mutually independent copies of the interarrival random times T and symbol * 
denotes the convolution operation 

t 

vp * vp(t) = J ip{t — r)ip(r)dT. 
o 

For the standard Poisson process, 

r n {t) = ^^r^e-^, (10) 

and according to the Central Limit Theorem 

¥ n) (t) = (^i/^* n (n/fi + t^i/fi) -^e-* 2/2 , n -)• oo. 

As numerical calculations show, ^"^(t) practically reaches its limit form already by n = 10 
(Fig. 2). 




-3-2-10 1 2 

t 



Figure 2: Rescaled arrival time distributions for the standard Poisson process 
{u = l, n = 1,2,3,5,10,30). 

In case of the fPp, 

oo 

ET = J ip u (t)tdt = oo 
o 
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and the Central limit theorem is not applicable. Applying the Generalized limit theorem (see, 
for example, [Uchaikin & Zolotarev, 1999]), we obtain: 



n 



l/u 



71 



o *n 

n l ' v i) v {tn l l v ) => g<r)(t), n ->• oo, 



where 



i> v (t) = Mt)l=i = t u - 1 E^ l/ (-f). 



Computing this multiple integrals can be performed by Monte Carlo technique. Taking 
fi = 1 and observing that ^™{t) is the probability density of the renormalized sum (T\ + 

o 

T-i + . . . + T n )/n < v of independent random variable, distributed according to ip u (t), we could 
directly simulate this sum by making use of the algorithm given by Eq. (jHJ) and construct the 
corresponding histogram. However, the left tail of the densities is too steep for this method, and 
we applied some modification of Monte Carlo method based on the partial analytical averaging 
of the last term. 




1 e-02 



1 e+02 



Figure 3: Rescaled arrival time distributions for fPp {y = 1/2; n = 1, 3, 10, and 30). 

By making use of this modification, we computed the distributions ^f^(t) for various n and 
v. An example of these results is represented in Fig. 3. 



6 The fractional Poisson distribution 

Now we consider another random variable: the number of events (pulses) N(t) arriving during 
the period t. According to the theory of renewal processes 

p n (t) = P(N(t) =n) = p(y / T j >t ] j-P (j^Tj >tj, n = 0, 1,2, . . . 

and the following system of integral equations for p n (t) takes place: 

oo t 

Pn(t) = S n0 J 7/j u (r)dT + [1 - 5 n0 ] J ip v (t - T)p n -i(r)dr, n = 0, 1, 2, . . . 

t 



After the Laplace transform with respect to time, we obtain 

X U p n (X) = -flp n (\) + yUPn-l(A) + y-^no, 

n = 0,1,2,..., p_i = 0. 
The inverse Laplace transform yields: 

D v tPn (t) = A/[p„-i(t) - p«(*)] + T( l_ u ^ no, < 1/ < 1. (11) 

This is a master equation system for the fractional Poisson processes. When v — > 1 it becomes 
the well known system for the standard Poisson process: 

= //[p n -i(t) - Pn (t)] + 5(t)5 n0 . (12) 
System ffTTl) produces for the generating function 

oo 

g(u,t) = J2n n Pn{t) (13) 

n=0 



the following equation: 



D^g{u, t) = fi(u - l)g(u, t) + ^ (14) 



When z/ — )■ 1 it becomes the well known equation for the standard Poisson process: 

dflr(ti,t) 



ix{u-l)g{u,t) + 5{t). (15) 



Comparing ( ITT]) with (112)1 and (JIl)) with (115)1 . one can observe that the equations for standard 
processes are generalized to the equations for correspondent fractional processes by means of 
replacement of the operator d/dt with Q D^ and of right side the term S(t) with t~ v /T(l — v). 
The solution to Eq. (j!4j) is of the form 

oo n 

Applying the binomial formula to each term of the sum and interchanging the summations, one 
can rewrite it as the series 



n=0 



a n ~ (m + n)!(-o) m 
n! n m\r(v(mk + n) + 1) 

m=0 v v ' y . 



(16) 



Comparing (fTBj) with (|T5|) yields 



a n ^ (m + n)! (-a) m 
PnW ~ nF m! r((m + n)i/+ 1)' 

This distribution, which becomes the Poisson one when v = 1 can be considered as its frac- 
tional generalization, called fractional Poisson distribution. The correspondent mean value and 
variance are given by 

fit u 



W)) 



~V(u + l) 
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and 

a\t) = (N(t)){l + (N(t))[2^uB(u, 1/2) - 1]}, 

where 

1 

B(ai,a 2 ) = Jx ai -\l - xY^dx 
o 

is the beta-function. 

7 Limit fractional Poisson distributions 

In case of the standard Poisson process, the probability distribution for random number N(t) 
of events follows the Poisson law with (N(t)) = fit = n which approaches to the normal one 
at large n. Introducing normalized random variable Z = N(t)/n and quasicontinuous variable 
z = n/n, one can express the last fact as follows: 



n" 2 



f ^ n)=n w^f n ~ 



n f (z-lf 
exp 




2tt L { 2/n J 

as n — > oo. In the limit case n — > oo the distribution of Z becomes degenerated one: 

Jim f(z;n) = S(z — 1). 

n— >oo 

Considering the case of fPp, we pass from the generating function to the Laplace characteristic 
function 

g(u,t) = E v (fit v (u-l)) = E v (nr(v + l)(u-l)). 
Introducing a new parameter A = — nlnu we get 

Eu N ® = Ee~ xz = E u (nT(u + l)(e" A/ " - 1)). 

At large n relating to large time t, 



Ee- xz = J e- Xz U{z)dz~E v {-}!), 
o 



A' = XT(u+l). 

Comparison of this equation with formula (6.9.8) of the book [Uchaikin & Zolotarev, 1999] 

E (-AO = Jj^l s ^ (x -u ndx = Wr+Mt M ( ) dz 

v J x w/»s {x > ax J e vz wi» 9 [n,, + D1-V- 







shows that the random variable Z has the non-degenerated limit distribution at t — > oo (see 
also [Uchaikin, 1999]): 

f (z- n) -+f(z)= [r( " + 1)]V % (-) ( Z "' V ) (17) 
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with moments 



feN [T(l + u)] k V(l + k) 



r(i + kv) 



By making use of series for g( v \ we obtain 

oo {~z) k 

UZ) = to k\T(l-(k+l)u)[T(u+l)]^- 



When z — >■ 0, 



/„(*) -> /„(()) 



1 



r(i + i/)r(i - 1/) z/7T 

It is also worth to note, that (Z°) = 1, (Z 1 ) = 1 and (Z 2 ) = 2z/B(z/, 1 + z/), so that the limit 
relative fluctuations are given by 



5 V = a m /{N) = J2vB(v, l+u)-l 



In particular cases 



5 = 1, 6 X = 0, 5 1/2 = Jtt/2 - 1 




Figure 4: Limit distributions ((T7D for z/ = 0.1(0.1)0.9 and 0.95. 
For v = 1/2, one can obtain an explicit expression for f v {z) : 

f 1/2 (z) = -e- z2 '\ z > 0. 

7T 

The family of this limit distributions are plotted in Fig. 4. 



8 Fractional Furry process 

Let us pass to the branching processes and consider its simplest case, when each particle 
converts into two identical ones at the end of its waiting time, distributed with density il) v (t). 
The process begins with one particle at t = and the first arrival time has the same distribution 
density ip v {t). When v = 1, the process is called the Furry process (Fp), therefore, in case of 
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v < 1 we can call it the fractional Furry process (fFp). The following integral equations govern 
the fFp: 

°° * 71-1 

Pn(t) = Snl J ^u(r)dT + [1 - 5 n0 - 5 n l] / lj) v (t - t) ^ p k (t )p n -k (t ) G?T, 71 = 1, 2, . . . 
t k=1 

Following the same way as before, we obtain 

"71-1 

D v t p n (t) = /i ^2pk{t)p n -k(t) -Pn{t) + 
k=l 



ri-i/ 



-<*nl, < 1/ < 1. 



The solution of this equation in case of v = 1 is well known: it is represented by the geometrical 
distribution 



Pnit) = e-# 



1-e 



71-1 



n = 1,2,3, 




v=1.0 



0.01 




v=0.8 



0.01 




v=0.5 



0.01 




v=0.2 



0.01 



Figure 5: Monte Carlo calculation of f v {z) for t = 5 and z/ = 1.0, 0.8, 0.5, 0.2 (histograms) by 
comparison with hypothetical distribution ([TBI (smooth lines). 

As to fFp for v < 1, we did not manage to derive the corresponding distribution from the 
fractional equation in a closed analytical form. The reason of the trouble lies in nonlinearity of 
the equation in case of branching. The only characteristics, the mean number of particles at 
time t has been found and expressed through the Mittag-Leffler function: 

(N(t)) = E v (jjf). 

All other results have been obtained by means of Monte Carlo simulation using the algorithm 
described above. 
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Observe, that in contrast to the fPp, the limit distribution of the normalized random variable 
Z in case of fFp is not degenerated. In particular, for the standard Furry process 



f(z) = Jim npn Z ([i Inn) = e z . 

n—>oo 

One could to suppose that in fractional case the "standard exponential function" is replaced 
with its fractional analogue 

/„(*) = z»- l E v , v {-z v ). (18) 
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Figure 6: \ 2 Goodness-of fit Test. 

Direct comparison of Monte Carlo data with formula ( Fl8|) (Fig. 5) allows to propose that 
they coincide at large t, and the % 2 goodness of fit analysis confirms this hypothesis (Fig. 6). 



9 Concluding remarks 

Considering the fractional Poisson process as an example of integer-valued fractional processes, 
one can suppose that the use of a-stable densities may occur very useful both for theoretical 
investigations and numerical simulations. Another example of integer-valued fractional pro- 
cesses, Furry branching process, has been too considered. We are planning to continue this 
work by analyzing binomial, negative binomial and some other integer-valued processes which 
can be useful for description of stochastic phenomena in laser physics, quantum optics and even 
in quantum chromodynamics i.e. quark-gluon plasma statistics [Botet & Ploszajczak, 2002]. 

This work is supported by Russian Foundation for Basic Research (project 07-01-00517). 
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